function F = cal_model_moment(param, const, eqm)

    % Read-in
    sigma           = param.sigma;
    beta_v          = param.beta_v;
    beta_m          = param.beta_m;    
    d               = param.d;
    Q               = param.Q;    
    gamma           = const.gamma;
    phi_v           = const.phi_v;
    phi_y           = const.phi_y;
    lnzQ            = const.lnzQ;        
    density_omega   = const.density_omega;    
    C_v0            = const.C_v0; 
    C_m             = const.C_m;    
    lnwage_grid     = const.lnwage_grid;  
    q_indicator     = eqm.q_indicator;
    nq              = eqm.nq;
    nd              = eqm.nd;
    quintileflag    = eqm.quintileflag;   
    E_zx0           = eqm.integ.E_zx0;
    E_zx1           = eqm.integ.E_zx1;
    E_zv0           = eqm.integ.E_zv0;
    E_zv1           = eqm.integ.E_zv1;  
    E_zm0           = eqm.integ.E_zm0;        
    E_zm1           = eqm.integ.E_zm1;           
    ExportProbQ     = eqm.ExportProbQ;   
    C_v1            = eqm.C_v1;    
    Pi0             = eqm.Pi0;
    Pi1             = eqm.Pi1;
    V               = eqm.V;
    M               = eqm.M;
    theta_v         = eqm.theta_v;
    theta_m         = eqm.theta_m; 
    X_m             = eqm.X_m;
    Psig            = eqm.Psig;
    c               = eqm.c;    
    rv_ads          = eqm.rv_ads;
    rv_rev          = eqm.rv_rev;   
    
    
    
%% Model Moments: Exporter Share and Export Intensity

    % Fraction of exporters 
    exporter_share_q = density_omega'*(ExportProbQ.*q_indicator)./nq;
    exporter_share_q(nq==0)=0;
    exporter_share_d = (sum(repmat(exporter_share_q.*nq,d,1).*quintileflag,2)./sum(repmat(nq,d,1).*quintileflag,2))';
    exporter_share = sum(exporter_share_q.*nq);  
        
    % Export intensity of exporters
    export_intensity_q = 1-rv_rev;
    export_intensity_d = (sum(repmat((1-rv_rev).*exporter_share_q.*nq,d,1).*quintileflag,2)./sum(repmat(exporter_share_q.*nq,d,1).*quintileflag,2))';
    export_intensity = sum((1-rv_rev).*nq);
    
    % Useful notation
    nq_exp = nq.*exporter_share_q;
    nq_nonexp = nq.*(1-exporter_share_q);
    
    % Store
    model.exporter_share_q      = exporter_share_q;
    model.exporter_share_d      = exporter_share_d;
    model.exporter_share        = exporter_share;
    model.export_intensity_q    = export_intensity_q;
    model.export_intensity_d    = export_intensity_d;
    model.export_intensity      = export_intensity;
    
    
    
%% Model Moments: Mean In-degree and Out-degree

    % Useful notation
    Vbar = C_v0.*Pi0.^(1/beta_v).*E_zv0 + rv_ads.*C_v1.*Pi1.^(1/beta_v).*E_zv1;
        
    % Mean in-degree by q
    mean_indegree_q = M.*theta_m./nq;
    mean_indegree_q_exp = theta_m./nq_exp.*C_m.*Pi1.^(1/beta_m).*E_zm1;
    mean_indegree_q_nonexp = theta_m./nq_nonexp.*C_m.*Pi0.^(1/beta_m).*E_zm0;
    
    % Mean out-degree by q
    mean_outdegree_q = Vbar.*sum(phi_v.*repmat(theta_v',1,Q),1)./nq;   
    mean_outdegree_q_exp = rv_ads.*C_v1.*Pi1.^(1/beta_v).*E_zv1.*sum(phi_v.*repmat(theta_v',1,Q),1)./nq_exp;
    mean_outdegree_q_nonexp = C_v0.*Pi0.^(1/beta_v).*E_zv0.*sum(phi_v.*repmat(theta_v',1,Q),1)./nq_nonexp;
    
        % Adjustment
        mean_indegree_q(nq==0)=0;    
        mean_indegree_q_exp(nq==0)=0;    
        mean_indegree_q_nonexp(nq==0)=0;   
        mean_outdegree_q(nq==0)=0;     
        mean_outdegree_q_exp(nq==0)=0;   
        mean_outdegree_q_nonexp(nq==0)=0;

        % Check
        %check1 = mean(abs(mean_indegree_q_nonexp.*nq_nonexp + mean_indegree_q_exp.*nq_exp - mean_indegree_q.*nq));
        %check2 = mean(abs(mean_outdegree_q_nonexp.*nq_nonexp + mean_outdegree_q_exp.*nq_exp - mean_outdegree_q.*nq));

    % Mean in-degree by quintile
    mean_indegree_d =  (sum(repmat(mean_indegree_q.*nq,d,1).*quintileflag,2)./sum(repmat(nq,d,1).*quintileflag,2))';
    mean_indegree_d_exp = (sum(repmat(mean_indegree_q_exp.*nq_exp,d,1).*quintileflag,2)./sum(repmat(nq_exp,d,1).*quintileflag,2))';
    mean_indegree_d_nonexp = (sum(repmat(mean_indegree_q_nonexp.*nq_nonexp,d,1).*quintileflag,2)./sum(repmat(nq_nonexp,d,1).*quintileflag,2))';
    
    % Mean out-degrees by quintile
    mean_outdegree_d =(sum(repmat(mean_outdegree_q.*nq,d,1).*quintileflag,2)./sum(repmat(nq,d,1).*quintileflag,2))';  
    mean_outdegree_d_exp = (sum(repmat(mean_outdegree_q_exp.*nq_exp,d,1).*quintileflag,2)./sum(repmat(nq_exp,d,1).*quintileflag,2))';  
    mean_outdegree_d_nonexp = (sum(repmat(mean_outdegree_q_nonexp.*nq_nonexp,d,1).*quintileflag,2)./sum(repmat(nq_nonexp,d,1).*quintileflag,2))';  
         
    % Mean in-degrees for all
    mean_indegree = sum(mean_indegree_q.*nq);
    mean_indegree_exp = sum(mean_indegree_q_exp.*nq_exp)./sum(nq_exp);
    mean_indegree_nonexp = sum(mean_indegree_q_nonexp.*nq_nonexp)./sum(nq_nonexp);
   
    % Mean out-degree for all
    mean_outdegree = sum(mean_outdegree_q.*nq);
    mean_outdegree_exporter = sum(mean_outdegree_q_exp.*nq_exp)./sum(nq_exp);
    mean_outdegree_nonexp = sum(mean_outdegree_q_nonexp.*nq_nonexp)./sum(nq_nonexp);
       
    % Store
    model.mean_indegree_q           = mean_indegree_q;
    model.mean_indegree_q_exp       = mean_indegree_q_exp;
    model.mean_indegree_q_nonexp    = mean_indegree_q_nonexp;
    model.mean_indegree_d           = mean_indegree_d;
    model.mean_indegree_d_exp       = mean_indegree_d_exp;
    model.mean_indegree_d_nonexp    = mean_indegree_d_nonexp;     
    model.mean_indegree             = mean_indegree;  
    model.mean_indegree_exp         = mean_indegree_exp;
    model.mean_indegree_nonexp      = mean_indegree_nonexp;    
    model.mean_outdegree_q          = mean_outdegree_q;
    model.mean_outdegree_q_exp      = mean_outdegree_q_exp;
    model.mean_outdegree_q_nonexp   = mean_outdegree_q_nonexp;
    model.mean_outdegree_d          = mean_outdegree_d;
    model.mean_outdegree_d_exp      = mean_outdegree_d_exp;  
    model.mean_outdegree_d_nonexp   = mean_outdegree_d_nonexp;
    model.mean_outdegree            = mean_outdegree;
    model.mean_outdegree_exp        = mean_outdegree_exporter;
    model.mean_outdegere_nonexp     = mean_outdegree_nonexp;
    
    
 
%% Model Moments: Mean and Variance of Sales

    % Useful moments
    E_lnz0 = density_omega'*(lnzQ.*(1-ExportProbQ).*q_indicator)./nq;
    E_lnz1 = density_omega'*(lnzQ.*ExportProbQ.*q_indicator)./nq;
    E2_lnz0 = density_omega'*(lnzQ.^2.*(1-ExportProbQ).*q_indicator)./nq;
    E2_lnz1 = density_omega'*(lnzQ.^2.*ExportProbQ.*q_indicator)./nq;
    
        % Adjustment
        E_lnz0(nq==0)=0;
        E_lnz1(nq==0)=0;
        E2_lnz0(nq==0)=0;
        E2_lnz1(nq==0)=0;    
        
    % Mean of log sales by q
    mean_lnsales_q = (sigma-1)*gamma*E_lnz0 + log(Pi0).*(1-exporter_share_q)...
                   + (sigma-1)*gamma*E_lnz1 + log(Pi1).*exporter_share_q;                        
    mean_lnsales_q_exp = (sigma-1)*gamma*E_lnz1.*nq./nq_exp + log(Pi1);
    mean_lnsales_q_nonexp = (sigma-1)*gamma*E_lnz0.*nq./nq_nonexp + log(Pi0);          
      
        % Adjustment
        mean_lnsales_q(nq==0)=0;
        mean_lnsales_q_exp(nq==0)=0;
        mean_lnsales_q_nonexp(nq==0)=0;    
    
    % Mean of log sales by quintile           
    mean_lnsales_d = (sum(repmat(mean_lnsales_q.*nq,d,1).*quintileflag,2)./sum(repmat(nq,d,1).*quintileflag,2))';
    mean_lnsales_d_exp = (sum(repmat(mean_lnsales_q_exp.*nq_exp,d,1).*quintileflag,2)./sum(repmat(nq_exp,d,1).*quintileflag,2))';
    mean_lnsales_d_nonexp = (sum(repmat(mean_lnsales_q_nonexp.*nq_nonexp,d,1).*quintileflag,2)./sum(repmat(nq_nonexp,d,1).*quintileflag,2))';
    
    % Mean of log sales
    mean_lnsales = sum(mean_lnsales_q.*nq);
    mean_lnsales_exp = sum(mean_lnsales_q_exp.*nq_exp)./sum(nq_exp);
    mean_lnsales_nonexp = sum(mean_lnsales_q_nonexp.*nq_nonexp)./sum(nq_nonexp);
    
    % Variance of log sales by q
    var_lnsales_q = ((sigma-1)*gamma)^2*(E2_lnz0) + 2*(sigma-1)*gamma*log(Pi0).*E_lnz0 + log(Pi0).^2.*(1-exporter_share_q)...
                    +((sigma-1)*gamma)^2*(E2_lnz1) + 2*(sigma-1)*gamma*log(Pi1).*E_lnz1 + log(Pi1).^2.*exporter_share_q...
                   -(mean_lnsales_q).^2;
            
    % Variance of log sales by quintile           
    var_lnsales_d = (sum(repmat((var_lnsales_q+mean_lnsales_q.^2).*nq,d,1).*quintileflag,2)./sum(repmat(nq,d,1).*quintileflag,2))'-mean_lnsales_d.^2;   
        
    % Variance of log sales
    var_lnsales = sum((var_lnsales_q+mean_lnsales_q.^2).*nq)-mean_lnsales.^2; 
    
    
    % Store
    model.mean_lnsales_q = mean_lnsales_q;
    model.mean_lnsales_q_exp = mean_lnsales_q_exp;
    model.mean_lnsales_q_nonexp = mean_lnsales_q_nonexp;
    model.mean_lnsales_d = mean_lnsales_d;
    model.mean_lnsales_d_exp = mean_lnsales_d_exp;
    model.mean_lnsales_d_nonexp = mean_lnsales_d_nonexp;
    model.mean_lnsales = mean_lnsales;
    model.mean_lnsales_exp = mean_lnsales_exp;
    model.mean_lnsales_nonexp = mean_lnsales_nonexp;
    model.var_lnsales_q = var_lnsales_q;
    model.var_lnsales_d = var_lnsales_d;
    model.var_lnsales = var_lnsales;    
    


%% Model Moments: Network Sales Share
        
    % Mean network sales by q
    mean_network_sales_q = (Pi0.*E_zx0 + rv_rev.*Pi1.*E_zx1)./nq;
    mean_network_sales_q_exp = rv_rev.*Pi1.*E_zx1./nq_exp;
    mean_network_sales_q_nonexp = Pi0.*E_zx0./nq_nonexp; 
    
        % Adjustment
        mean_network_sales_q(nq==0)=0;
        mean_network_sales_q_exp(nq==0)=0;
        mean_network_sales_q_nonexp(nq==0)=0;
    
    % Mean network sales by quintile
    mean_network_sales_d = (sum(repmat(mean_network_sales_q.*nq,d,1).*quintileflag,2)./sum(repmat(nq,d,1).*quintileflag,2))'; 
    mean_network_sales_d_exp = (sum(repmat(mean_network_sales_q_exp.*nq_exp,d,1).*quintileflag,2)./sum(repmat(nq_exp,d,1).*quintileflag,2))'; 
    mean_network_sales_d_nonexp = (sum(repmat(mean_network_sales_q_nonexp.*nq_nonexp,d,1).*quintileflag,2)./sum(repmat(nq_nonexp,d,1).*quintileflag,2))';     

    % Network sales share by q
    network_sales_share_q = mean_network_sales_q.*nq/sum(mean_network_sales_q.*nq);
    network_sales_share_q_exp = mean_network_sales_q_exp.*nq_exp/sum(mean_network_sales_q_exp.*nq_exp);        
    network_sales_share_q_nonexp = mean_network_sales_q_nonexp.*nq_nonexp/sum(mean_network_sales_q_nonexp.*nq_nonexp);        
        
        % Adjustment
        network_sales_share_q(nq==0)=0;
        network_sales_share_q_exp(nq==0)=0;
        network_sales_share_q_nonexp(nq==0)=0; 
    
    % Network sales share by quintile
    network_sales_share_d = mean_network_sales_d.*nd/sum(mean_network_sales_d.*nd);        
    network_sales_share_d_exp = mean_network_sales_d_exp.*nd.*exporter_share_d/sum(mean_network_sales_d_exp.*nd.*exporter_share_d);        
    network_sales_share_d_nonexp = mean_network_sales_d_nonexp.*nd.*(1-exporter_share_d)/sum(mean_network_sales_d_nonexp.*nd.*(1-exporter_share_d));        
        
    % Store 
    model.network_sales_share_q         = network_sales_share_q;
    model.network_sales_share_q_exp     = network_sales_share_q_exp;
    model.network_sales_share_q_nonexp  = network_sales_share_q_nonexp;
    model.network_sales_share_d         = network_sales_share_d;    
    model.network_sales_share_d_exp     = network_sales_share_d_exp;
    model.network_sales_share_d_nonexp  = network_sales_share_d_nonexp;
        
    
    
%% Model Moments: Aggregate Network Moments

    % Supplier and customer shares
    common = mysum(repmat(theta_v',1,Q).*phi_v.*repmat(Vbar,Q,1), quintileflag, param);
    share_indegree_d = common./repmat(sum(common,2),1,d);
    share_outdegree_d = common./repmat(sum(common,1),d,1);

    % Expenditure and sales share
    common = mysum(repmat((theta_m./V.*c.^(sigma-1).*X_m)',1,Q).*phi_y.*phi_v.*repmat(Psig,Q,1), quintileflag, param);
    share_spend_d = common./repmat(sum(common,2),1,d);
    share_sales_d = common./repmat(sum(common,1),d,1);
    
    % Store
    model.share_indegree_d  = share_indegree_d;
    model.share_outdegree_d = share_outdegree_d;
    model.share_spend_d     = share_spend_d;
    model.share_sales_d     = share_sales_d;
    
        
    
%% Model Moments: Average Log Wage of Suppliers

    % Average supplier wage by q (same for all/exporter/nonexporter)    
    avg_unwgt_lnwageS_q = 1./V.*sum(phi_v.*repmat((log(param.w)+lnwage_grid).*Vbar, Q,1), 2)';    
    avg_wgt_lnwageS_q = theta_m.*c.^(sigma-1)./V.*sum(phi_y.*phi_v.*repmat((log(param.w)+lnwage_grid).*Psig, Q,1),2)';    
    
    % Average supplier wage by quintile
    avg_wgt_lnwageS_d = (sum(repmat(avg_wgt_lnwageS_q.*nq, d, 1).*quintileflag, 2)./sum(repmat(nq,d,1).*quintileflag,2))';  
    avg_wgt_lnwageS_d_exp = (sum(repmat(avg_wgt_lnwageS_q.*nq_exp, d, 1).*quintileflag, 2)./sum(repmat(nq_exp,d,1).*quintileflag,2))';  
    avg_wgt_lnwageS_d_nonexp = (sum(repmat(avg_wgt_lnwageS_q.*nq_nonexp, d, 1).*quintileflag, 2)./sum(repmat(nq_nonexp,d,1).*quintileflag,2))';  
    avg_unwgt_lnwageS_d = (sum(repmat(avg_unwgt_lnwageS_q.*nq, d, 1).*quintileflag, 2)./sum(repmat(nq,d,1).*quintileflag,2))';   
    avg_unwgt_lnwageS_d_exp = (sum(repmat(avg_unwgt_lnwageS_q.*nq_exp, d, 1).*quintileflag, 2)./sum(repmat(nq_exp,d,1).*quintileflag,2))';   
    avg_unwgt_lnwageS_d_nonexp = (sum(repmat(avg_unwgt_lnwageS_q.*nq_nonexp, d, 1).*quintileflag, 2)./sum(repmat(nq_nonexp,d,1).*quintileflag,2))';   
    
    % Average supplier wage
    avg_wgt_lnwageS = sum(avg_wgt_lnwageS_q.*nq);
    avg_wgt_lnwageS_exp = sum(avg_wgt_lnwageS_q.*nq_exp)./sum(nq_exp);
    avg_wgt_lnwageS_nonexp = sum(avg_wgt_lnwageS_q.*nq_nonexp)./sum(nq_nonexp);
    avg_unwgt_lnwageS = sum(avg_unwgt_lnwageS_q.*nq);
    avg_unwgt_lnwageS_exp = sum(avg_unwgt_lnwageS_q.*nq_exp)./sum(nq_exp);
    avg_unwgt_lnwageS_nonexp = sum(avg_unwgt_lnwageS_q.*nq_nonexp)./sum(nq_nonexp);
    
    % Store    
    model.avg_wgt_lnwageS_q           = avg_wgt_lnwageS_q;   
    model.avg_wgt_lnwageS_d           = avg_wgt_lnwageS_d;
    model.avg_wgt_lnwageS_d_exp       = avg_wgt_lnwageS_d_exp;
    model.avg_wgt_lnwageS_d_nonexp    = avg_wgt_lnwageS_d_nonexp;
    model.avg_wgt_lnwageS             = avg_wgt_lnwageS;
    model.avg_wgt_lnwageS_exp         = avg_wgt_lnwageS_exp;
    model.avg_wgt_lnwageS_nonexp      = avg_wgt_lnwageS_nonexp;   
    model.avg_unwgt_lnwageS_q         = avg_unwgt_lnwageS_q;
    model.avg_unwgt_lnwageS_d         = avg_unwgt_lnwageS_d;
    model.avg_unwgt_lnwageS_d_exp     = avg_unwgt_lnwageS_d_exp;
    model.avg_unwgt_lnwageS_d_nonexp  = avg_unwgt_lnwageS_d_nonexp;
    model.avg_unwgt_lnwageS           = avg_unwgt_lnwageS;
    model.avg_unwgt_lnwageS_exp       = avg_unwgt_lnwageS_exp; 
    model.avg_unwgt_lnwageS_nonexp    = avg_unwgt_lnwageS_nonexp;  
    model.avg_unwgt_lnwageS_d4    = model.avg_unwgt_lnwageS_d(2:5) - model.avg_unwgt_lnwageS_d(1);
    model.avg_wgt_lnwageS_d4      = model.avg_wgt_lnwageS_d(2:5) - model.avg_wgt_lnwageS_d(1);    
    
    
    
%% Model Moments: Wage Regression Coefficients

    % Decompose weighted average log wage of suppliers
    avg_wgt_lnwageS_ext_q = avg_unwgt_lnwageS_q;
    avg_wgt_lnwageS_int_q = avg_wgt_lnwageS_q - avg_unwgt_lnwageS_q;    
    
    % Wage regression
    reg_W = diag(nq);
    reg_x = [ones(Q,1), param.w.*lnwage_grid'];
    reg_y_ext = avg_wgt_lnwageS_ext_q';
    reg_y_int = avg_wgt_lnwageS_int_q';
    reg_y_sum = avg_wgt_lnwageS_q';
    beta_ext = (reg_x'*reg_W*reg_x)\(reg_x'*reg_W*reg_y_ext);
    beta_int = (reg_x'*reg_W*reg_x)\(reg_x'*reg_W*reg_y_int);
    beta_sum = (reg_x'*reg_W*reg_x)\(reg_x'*reg_W*reg_y_sum);
    slope_ext = beta_ext(2);
    slope_int = beta_int(2);
    slope_sum = beta_sum(2);   

    % Store    
    model.avg_wgt_lnwageS_ext_q   = avg_wgt_lnwageS_ext_q;
    model.avg_wgt_lnwageS_int_q   = avg_wgt_lnwageS_int_q;
    model.slope_ext             = slope_ext;
    model.slope_int             = slope_int;
    model.slope_sum             = slope_sum;    
    
        
    
%% Match Quality Composition

    % Fraction of suppliers of higher (or same) quality
    flag_S_higher_q = zeros(param.Q, param.Q); %row is buyer, column is seller
    for i = 1:param.Q
        flag_S_higher_q(i,:) = [zeros(1,i-1), ones(1,param.Q-i+1)];
    end
    model.frac_S_higher_q = sum(flag_S_higher_q.*const.phi_v.*repmat(eqm.Vbar,param.Q,1),2)'./eqm.V;    
    
    
    % Fraction of customers of higher (or same) quality
    flag_B_higher_q = zeros(param.Q, param.Q); %row is buyer, column is seller
    for i = 1:param.Q
        flag_B_higher_q(:,i) = [zeros(i-1,1); ones(param.Q-i+1,1)];
    end
    tempmatrix = const.phi_v.*repmat(eqm.theta_v',1,param.Q);
    model.frac_B_higher_q = sum(flag_B_higher_q.*tempmatrix,1)./sum(tempmatrix,1);        
    
    
    
%% Number of Workers

    % Total wage bill by q
    model.total_wage_bill_q = (1-param.alpha_m-param.alpha_s)*(param.sigma-1)/param.sigma*(eqm.X0+eqm.X1);
    model.total_wage_bill_q_exp = (1-param.alpha_m-param.alpha_s)*(param.sigma-1)/param.sigma*eqm.X1;
    model.total_wage_bill_q_nonexp = (1-param.alpha_m-param.alpha_s)*(param.sigma-1)/param.sigma*eqm.X0;
        
    % Total number of workers by q
    model.total_num_worker_q = model.total_wage_bill_q./const.wage_grid;
    model.total_num_worker_q_exp = model.total_wage_bill_q_exp./const.wage_grid;
    model.total_num_worker_q_nonexp = model.total_wage_bill_q_nonexp./const.wage_grid;
    
        % Adjustment
        model.total_num_worker_q(nq==0) = 0;
        model.total_num_worker_q_exp(nq==0) = 0;
        model.total_num_worker_q_nonexp(nq==0) = 0;
        
        
    
%% RETURN
    
    F = model;

end